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ABSTRACT 

We present a novel, fast method to recover the density field through the statistics of the trans- 
mitted flux in high redshift quasar absorption spectra. The proposed technique requires the 
computation of the probability distribution function of the transmitted flux (Pp) in the Lya 
forest region and, as a sole assumption, the knowledge of the probability distribution function 
of the matter density field (Pa)- We show that the probability density conservation of the flux 
and matter density unveils a flux-density (P — A) relation which can be used to invert the 
Lya forest without any assumption on the physical properties of the intergalactic medium. 
We test our inversion method at z = 3 through the following steps: [i] simulation of a sample 
of synthetic spectra for which Pa, is known; [ii] computation of Pp; [iii] inversion of the Lya 
forest through the P — A relation. Our technique, when applied to only 10 observed spec- 
tra characterized by a signal-to noise ratio S / N > 100 provides an exquisite (relative error 
€a ;$ 12% in > 50% of the pixels) reconstruction of the density field in > 90% of the line of 
sight. We finally discuss strengths and limitations of the method. 
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1 INTRODUCTION 



The ultraviolet radiation emitted by a quasar can suffer resonant 
Lya scattering as it propagates through the intergalactic neutral hy- 
drogen. In this process, photons are removed from the line of sight 
resulting in an attenuation of the source flux, the so-called Gunn- 
Peterson effect (Gunn & Peterson 1965). At z ~ 3, the neutral hy- 
drogen number density is sufficiently small to allow to resolve the 
wide series of absorption lines, which give rise to the so-called Lya 
forest in quasar absorption spectra. Initially, the absorption lines 
observed in quasar spectra were believed to be produced by an in- 
tergalactic population of pressure-confined clouds, photoionized by 
the integrated quasar flux (Sargent et al. 1980). Nowaday absorbers 
of the Lya forest are generally associated to large-scale neutral hy- 
drogen density fluctuations in the warm photoionized IGM. 
This scenario, firstly tested through analytical calculations (Bond 
et al. 1988; Bi et al. 1992), has received an increased consensus 
thanks to the detection of a clustering signal both along individual 
(Cristiani et al. 1995; Lu et al. 1996; Cristiani et al. 1997; Kim et al. 
2001) and close pairs of quasar lines of sight, resulting in an esti- 
mated size for the absorbers which varies between few hundreds of 
kpc (Smette et al. 1992; D'Odorico et al. 1998; Rauch et al. 2001; 
Becker et al. 2004) to few Mpc (Rollinde et al. 2003; Coppolani 
et al. 2006; D'Odorico et al. 2006). Hydrodynamical cosmological 
simulations lend further support to this picture, having established 
the existence of a tight connection between the Hi density field and 
the dark matter distribution, at least on scales larger than the leans 



lenght of the IGM (Cen et al. 1994; Petitjean, et al. 1995; Zhang 
et al. 1995; Hernquist et al. 1996; Miralda-Escude' et al. 1996). As 
a consequence of this result, the Lya forest has been proposed as 
a powerful method to study the clustering properties of the dark 
matter density field and to measure its power spectrum (Croft et 
al. 1998; Nusser & Haehnelt 1999; Pichon et al. 2001; Rollinde et 
al. 2001; Zhan et al. 2003; Viel et al. 2004; McDonald et al. 2005; 
Zaroubi et al. 2006; Saitta et al. 2008; Bird et al. 2010). 
In this paper, we focus our attention on the reconstruction of the 
density field along the line of sight (LOS), which represents the 
starting point for the dark matter power spectrum measurement. 
We present a novel, fast method to recover the density field along 
the LOS towards z ~ 3 quasar which can be considered alternative 
and/or complementary to the existing techniques in the literature 
mentioned above which will be further discussed in this work. 



2 MODELLING THE FOREST 

Let us discretize the LOS towards a distant quasar in a number of 
pixels iVpi x . The transmitted flux due to Lya absorption in the IGM 
at a given pixel i is computed from the usual relation 



F(i) 



-r(i) 



(1) 



where r(i) is the Lya absorption optical depth: 
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Figure 1. Left panel: Probability distribution function of the transmitted flux obtained from 10 synthetic quasar absorption spectra at z = 3. Right panel: 
Probability distribution function of the density field at z = 3, as predicted by the lognormal model. The thick line shows P/± of the recovered density field. The 
probability of finding a pixel characterized by a transmitted flux above a given F maJ( (cyan shaded region on the left-hand panel) is equal to the probability of 
finding an overdensity below Ab (cyan shaded region on the right-hand panel). Analogously, the probability of finding a pixel characterized by a transmitted 
flux below a given F nl - ln (gray shaded region on the left-hand panel) is equal to the probability of finding an overdensity above Ad (gray shaded region on 
the right-hand panel). In particular, at each pixel characterized by a transmitted flux F» it can be associated an overdensity A* such that eq. (8) is satisfied. In 
this case, the probability of finding a pixel characterized by a transmitted flux F mm < F* < F max (red hatched region on the left-hand panel) is equal to the 
probability of finding an overdensity Ab < A* < A<j (red hatched region on the right-hand panel). 



-(i) = do 



Ax 



l + z(i) 



^2 n m (j)& a [vH.(i) - v-n(j)] ■ 



(2) 



3=1 



where c is the light speed, I a is the Lya cross section, Aa; is the 
comoving pixel size, z(i) is the redshift of pixel i, $ a is the Voigt 
profile for the Lya transition, vh is the Hubble velocity 1 , and nm 
is the neutral hydrogen fraction. 

The computation of nm generally assumes that the low density gas 
which gives rise to the Lya forest is approximately in local equi- 
librium between photoionization and recombination. By neglecting 
the presence of helium, and assuming that the IGM is highly ion- 
ized at the redshift of interest by a uniform ultraviolet background, 
the photoionization equilibrium condition provides us with the fol- 
lowing relation: 



nm(J) 



<*[T(j)] 



Wl + z(j)] 3 A(j)} 2 cxA0f 



(3) 



with 1.5 < P < 2.0 (Hui & Gnedin 1997). In eq. (3), Y is 
the photoionization rate of neutral hydrogen at a given redshift 
z(j), A(j) is the overdensity at the pixel j, no is the mean 
baryon number density at z = 0, and a[T(j)] oc T(j)~ ' 7 is 
the temperature-dependent radiative recombination rate. For quasi- 
linear IGM, where non-linear effects like shock-heating can be 
neglected, the temperature can be related to the baryonic density 
through a power-law relation (Hui & Gnedin 1997): 



T(j) = T A(j) 



7-1 



(4) 



where To is the IGM temperature at the mean density at a given 
redshift z(j) which depends on the reionization history of the Uni- 
verse, such as the slope 7 of the equation of state, and T. 
Besides the nm computation, the thermal state of the IGM affects 
the profile of the absorption lines. In fact, the Voigt profile which 
enters in eq. (2) is the convolution of a Gaussian and a Lorentzian 



The velocity of a Lya forest absorber is generally given by the sum of 
the Hubble and the peculiar velocities, i.e. v(i) = «h(0 + v pec (i)- In 
this work, we neglect peculiar velocities and we discuss this issue in the 
Discussion section. 



profile, since it takes into account the thermal and the natural broad- 
ening of the absorption line, respectively. For low column density 
regions, as the ones corresponding to Lya forest, thermal broaden- 
ing dominates, and the Voigt function reduces to a simple Gaussian 



$a[VH(i)— UhO')] 



exp 



vn(i) — vn(j) 



b(j) 



.(5) 



whose width is set by the temperature-dependent Doppler parame- 
ter b(j) — ^/2kBT(j)/m p , with fee and ra p being the Boltzmann 
constant and the proton mass, respectively. 

From the system of equations (l)-(5) the complex relation between 
the transmitted flux F in quasar absorption spectra and the under- 
lying density field A becomes evident. The determination of the 
F — A relation not only depends on the parameters describing the 
physical state of the IGM (r; To; 7), i.e. on the cosmic reionization 
history, but also on the convolution of the resulting nm with the 
profile of the absorption lines, which in turn is linked to the ther- 
mal state of the gas through the Doppler parameter b. In the next 
Section, we describe how the F — A relation can be efficiently in- 
ferred from a statistical analysis of the transmitted flux in quasar 
absorption spectra. 



3 FLUX-DENSITY RELATION 

Let us assume that the signal-to-noise ratio characterizing an ob- 
served quasar absorption spectrum is such that the maximum flux 
(minimum optical depth) which can be distinguished from full 
transmission is given by F max (r m i n ), and that the minimum de- 
tectable flux (maximum optical depth) is given by F m i n (r max ). By 
combining eq. (2)-(3)-(4), it results that r = t(z, ft, V, T , 7, A). 
This means that, at a given redshift z, for any /?, and for any given 
IGM thermal and ionization history, two characteristic overden- 
sities do exist: [i] Ab, where the subscript "b" means "bright", 
such that each Lya absorber sitting on A < Ab would provide 
r(Ab) < T m i n , i.e. full transmission of the quasar continuum 
(F > Fmax); [ii] Ad, where the subscript "d" means "dark", such 
that each overdensity A > Ad would provide r(Ad) > T max , 
therefore completely depleting the flux emitted by the source at the 
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Figure 2. Flux-density relation resulting from flux/density probability conservation is shown for different parameters defining the IGM physical properties, as 
labelled in the figures. The black squares represent results obtained from 10 synthetic spectra, while coloured lines denote the flux-density relations resulting 
by solving eq. (8). In the rightmost panel, the flux-density relation is plotted for all the set of parameters considered, which provide the same mean flux. 



location of the absorber (F < F m - m ). 

One of the statistics adopted to analyze quasar spectra is the prob- 
ability distribution function of the transmitted flux (Pf), which 
quantify the number of pixels characterized by a transmitted flux 
between F and F + dF. Let us assume that at a given overden- 
sity A* it is associated a unique value of the transmitted flux F*. 
In this case, Ab and Ad can be defined through the probability of 
finding a pixel characterized by a transmitted flux F > F max , and 
F < F m i n , respectively. Once Pf is computed from observed spec- 
tra, Ab and Ad can be determined from the following equations: 



i Mb 

P F dF = / P A dA , 

-Fmax J 

lin r+oo 

Pp dF — I P A dA . 



(6) 



(7) 



where Pa is the probability distribution function of the density 
field. Analogously to eq. (6) and eq. (7), to each pixel character- 
ized by a transmitted flux P m in < P* < Pmax an overdensity 
Ab < A* < Ad can be associated such that: 



Pp dF = / Pa dA , 

F. JA b 

or, equivalently, 



F, fAd 

PpdF = / P A dA . 



(8) 



(9) 



Note that eq. (6)-(7)-(8)-(9) are simply obtained by imposing the 
conservation of the flux/density probability densities. The inference 
of Ab, Ad, A* only depends on the assumed Pa and not on any 
assumption concerning thermal and ionization histories which are 
however encoded into the observed Pp, as we show in the next Sec- 
tion. By solving eq. (8), for each observed value of P„ in terms of 
A*, the flux-density relation of the Lya forest can be readily de- 
rived. Such F — A relation can be used to invert the Lya forest, 
therefore allowing to reconstruct the density field along the LOS. 



4 INVERSION PROCEDURE 

We simulate quasar absorption spectra at z = 3 through the model 
described in the previous section. In particular, for what concerns 
the baryonic density distribution along the LOS entering in eq. (3), 
we adopt the method described by Gallerani et al. (2006), whose 
main features are summarized as follows. The spatial distribution 



of the baryonic density field and its correlation with the peculiar 
velocity field are taken into account adopting the formalism intro- 
duced by Bi & Davidsen (1997). To enter the mildly non-linear 
regime which characterizes the Lya forest absorbers we use a log- 
normal model introduced by Coles & Jones (1991). 
The thermal and ionization state of the IGM are fixed in such 
a way that the mean transmitted flux simulated in a sample of 
lines of sight is F ~ 0.72, as obtained from real data (Songaila 
2004). We adopt the following set of parameters: [i] (r, To, 7) = 
(3.0, 1.5, 1.2); [ii] (3.2, 1.5,0.6); [iii] (4.1, 1.0, 1.2), where V is 
reported in units of 10 -12 s _1 , and To in units of 10 4 K. 
We add observational artifacts to the synthetic spectra: first, we 
smooth each simulated spectrum to a resolution R=36000, com- 
parable to the one of high-resolution spectrographs (e.g. HIRES); 
then we add noise to get a signal-to-noise ratio S/N=100; finally 
we rebin the noisy synthetic spectra to R= 10000, in order to re- 
move small flux fluctuations. Starting from a sample of 10 LOS, 
we compute the resulting Pf for the cases [i], [ii], [iii]. In Fig. 1, 
Pf for the case [i] is shown in the left panel, through a solid line. 
We assume the lognormal Pa which has been adopted for the syn- 
thetic spectra simulations, and which is described by a Gaussian 
with mean = —0.28, and sigma = 0.48 (Fig. 1, right panel, solid 
line). The choice of a lognormal Pa is only made for testing pur- 
poses, but when dealing with real data any other, perhaps more 
realistic, prescription for Pa can be implemented. For example, 
Miralda-Escude' et al. 2000 suggested a different functional form, 
whose agreement with data has not yet fully established (see dis- 
cussion in Becker et al, 2006). 

Starting from the computed Pf and the assumed Pa, we solve 
eq. (8), to establish the flux-density relation for each of the set of 
parameters considered. The resulting F — A relations 2 are shown 
through coloured lines in Fig. 2 for the case [i] (left-most panel, 
red line), [ii] (middle panel, green line), [iii] (right-most panel, blue 
line). In this figure, the filled squares and the solid black lines rep- 
resent the mean flux and the relative la error, respectively, corre- 
sponding to a given overdensity, as obtained from our simulations. 
It can be seen that the resulting F — A relations obtained by solv- 
ing eq. (8), agree very well with the relation between the mean flux 
and the overdensities predicted by the simulations. However, our 
method does not manage to take into account the dispersion of the 
F — A relations. 



2 In this work we adopt the cosmological parameters provided by Komatsu 
et al. (2010). Uncertainties on the cosmological parameters have negligible 



effect on the flux-density relation. 
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The reason for this is due to our assumption of a one-to-one corre- 
spondence between the transmitted flux observed at a given pixel 
i and the underlying matter density field at the corresponding red- 
shift z(i). Such an assumption would be valid if the profile of the 
absorption lines would be described by a Dirac Delta function (i.e. 
in the ideal case of an infinitely small natural broadening of the 
line), and in the case of a cold IGM (i.e. in the ideal case of no 
thermal broadening of the line). In this ideal case the cross section 
of the Lya transition is different from zero only in correspondence 
of the pixel i for which we compute the optical depth. In the real 
case, to the optical depth of the pixel i do contribute also the ad- 
jacent pixels. In particular, the higher is the IGM temperature (i.e. 
the larger is the doppler parameter), the larger is the number of ad- 
jacent pixels responsible for the optical depth of the pixel i. There- 
fore, an observed value of the transmitted flux F* at the pixel i can 
be produced by an isolated overdensity A > A*, as well as by 
a clustering of overdensities A < A*. Our method does not take 
into account the broadening of the absorption lines, therefore it can 
not provide the dispersion in the F — A relations; in the discussion 
section, we will analyze the effect of this limit of our method on 
the density field reconstruction. 

In Fig. 3, we show an example of the simulated density field for 
a random LOS (bottom panel, black line) and the corresponding 
transmitted fluxes (top panel) obtained for for the case [i] (red line), 
[ii] (green line), [iii] (blue line). In the bottom panel of Fig. 3 the 
originally adopted density field (black line) is compared with the 
recovered ones, as obtained by applying the F — A relation to 
the transmitted flux computed from the set of parameters [i] (red 
line), [ii] (green line), and [iii] (blue line). This figure shows that 
the method works quite well 3 for all the set of parameters consid- 
ered, i.e. almost independently from the physical properties of the 
IGM. This comes out by the fact that the F — A relations obtained 
by solving eq. (8) in terms of A* (coloured lines in Fig. 2) are 
not completely independent on the reionization history, i.e. differ- 
ent set of parameters affects differently the P&, therefore providing 
slightly different F — A relations. This is the reason why by in- 
verting the Lya forest we can associate the same density field to 
different values of the transmitted flux, obtaining a reconstruction 
method which is free of any assumption on the thermal and ioniza- 
tion state of the IGM. In the next section we discuss the strengths 
and the limits of our method in further details. 



5 DISCUSSION 

From the bottom panel of Fig. 3, it is evident that the method is 
overall very successful. Nevertheless some minor discrepancie be- 
tween the original and reconstructed density field remain. In partic- 
ular, the reconstructed density field is (1) smoother than the origi- 
nal one, and (2) it overestimates (underestimates) the original one 
in the most underdense (overdense) regions. 
The accuracy loss (1) arises from the adopted assumption of a one- 
to-one correspondence between the transmitted flux and the matter 
density. As already discussed in the previous section, this assump- 
tion does not allow us to take into account the clustering of the 
density field locally, or equivalently the thermal broadening of the 
absorption lines. This results in a reconstructed density field which 
is smoother than the original one. 

3 By computing the relative error ca of the recovered density field along 
> 10 lines of sight, we find that > 50% of the pixels are characterized by 
e A < 0.12. 
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Figure 3. Top panel: Example of synthetic quasar absorption spectra 
at z = 3 as obtained by adopt the following set of parameters: [i] 
(T,T ,7) = (3.0,1.5,1.2) (red line); [ii] (3.2,1.5,0.6) (green line); [iii] 
(4.1, 1.0, 1.2) (blue line), where Y is reported in units of 10~ 12 s _1 , and 
To in units of 10 4 K. Bottom panel: The dotted black line represents the 
distribution of the overdensities which produce the transmitted flux shown 
in the top panel. The red line (green/blue) show the recovered density field, 
as obtained by inverting the flux density relation corresponding to the set 
of parameters [i] [ii/iii] . The magenta line show the case of a high quality 
(S/N=1000) spectrum. 

As for limit (2), this is due to the fact that we are considering both a 
maximum flux (T max ) which can be distinguished from full trans- 
mission, and a minimum detectable flux (T m i n ). As explained in 
Sec. 3, these maximum and minimum transmitted fluxes translate 
into a minimum (Ab) and maximum (Ad) overdensities which can 
be recovered. Stated differently, in correspondence of the most un- 
derdense regions (A < Ab) we can only put an upper limit on 
the recovered density field, i.e. (A = Ab); analogously, for the 
most overdense regions (A > Ad) only a lower limit ocan be set, 
i.e. (A = Ad). For a fixed mean transmitted flux, the maximum 
and minimum recoverable overdensities depend on IGM proper- 
ties, and are more sensitive 4 to changes in 7 than to To: [i/ii] 
(A b ;A d ) = (0.09; 4.90); [iii] (A b ;A d ) = (0.11; 3.90). Typi- 
cal uncertainties on T mcan at z — 3 (< 10%, Songaila 2004) im- 
plies small variations on Ab and Ad. For example, in the case [i] 
a mean transmitted flux 10% higher [lower] than 0.72 results in 
(A b ; A d ) = (0.11; 6.10) [(A b ; A d ) = (0.08; 3.82)]. 
The existence of Ab and Ad allows us to recover the density field 
along > 92% of the simulated lines of sight. This result can be fur- 
ther improved by applying our method to very high quality spectra. 
Up to now we have assumed a S/N = 100, which means that we 
are sensitive at la level to flux variations up to the second decimal 
digit. Therefore, we have considered T max = 0.99 and T m i n = 
0.01. However, if we consider higher quality data (e.g. S/N=1000), 
the original density field is better recovered, since the completeness 
of the signal rises to > 98%, being (A b ; A d ) = (0.03; 6.60). In 
this case, we obtain a reconstructed signal characterized by a rela- 



4 For an inverted equation of state [iii], Ab is higher and Ad is lower than 
in the case of a "regular" one. In fact, as can be seen from Fig. 3, in the case 
[iii] (green line), the most underdense regions (A < Ab) tend to be more 
transmitting (F > T max ) and the most opaque regions ( A > Ad) to be 
even more absorbed (F < T m i n ) 
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five error e A < 12% in 51% of the pixels, 12% < ea < 42% in 
39% of the pixels, and 42% < e A < 100% in 10% of the pixels. 
We finally compare our method with previously developed tech- 
niques. One of the approaches generally adopted in the density field 
reconstruction is to define a physical model relating the observed 
flux to the underlying density field and then invert this relation. As 
shown in Sec. 2, any model which allows to predict the transmit- 
ted flux resulting from a given density field distribution, at a given 
redshift z, depends on several uncertain parameters (/?, Y, To 7). 
This approach has been pioneered by Nusser & Haehnelt (1999), 
hereafter NH99. These authors have developed a technique based 
on a model of the Lya forest, which allows to reconstruct the den- 
sity field along quasar sightlines through an iterative method. Such 
technique allows to correctly take into account the thermal broad- 
ening of the absorption lines, assuming that both To and 7, which 
affect the Doppler parameter, are known. 

The NH99 inversion method has been adopted by several authors 
to study the properties of the IGM and to measure the matter power 
spectrum from the Lya forest (e.g. Pichon et al. 2001; Rollinde 
et al. 2001; Zaroubi et al. 2006). However, it has been recognized 
that uncertainties on j3, 7, and To result in spurious biases in the 
recovered density field (NH99; Pichon et al. 2001; Rollinde et al. 
2001). The strong dependence of the matter power spectrum on 
the IGM properties, as inferred from Lya forest data, has been 
also confirmed by Bird et al. (2010). Saitta et al. (2008) have pro- 
posed another possible approach to reconstruct the density field 
through the transmitted flux in quasar absorption spectra, named 
FLO (i.e. From Lines to Overdensities). Also in this case the pa- 
rameters defining the IGM equation of state need to be assumed. 
Moreover, while the FLO method seems promising for the recon- 
struction of the most overdense regions (up to A ~ 30), underden- 
sities (A < 1) tend to be underestimated. 

Summarizing, the methods available so far in the literature tried to 
invert the F — A relation for which a large number of assumptions 
have to be made on the physical state of the IGM. We have shown 
that the F — A relation can be efficiently inferred from a statistical 
analysis of the transmitted flux in quasar absorption spectra dras- 
tically reducing the number of assumptions. As far as this work is 
concerned, we have not considered peculiar velocity effects which 
influence the quality of the recovery (NH99). However, we note 
that the same method described in NH99 for the reconstruction in 
real space can be applied to our recovered density field iteratively. 
We leave this to a future work. We finally note that so far we have 
been referring to quasar absorption spectra. However, the method 
proposed can be applied to GRB absorption spectra as well. 



6 CONCLUSIONS 

A novel, fast method to recover the density field through the statis- 
tics of the transmitted flux in high redshift quasar absorption spec- 
tra has been introduced. The proposed technique requires the com- 
putation of the probability distribution function of the transmit- 
ted flux (Pf) in the Lya forest region and, as a sole assumption, 
the knowledge of the probability distribution function of the mat- 
ter density field (Pa). We have shown that the probability density 
conservation of the flux and matter density unveils a flux-density 
(F — A) relation which can be used to invert the Lya forest, with- 
out any assumption on the properties of the IGM. 
Our inversion method has been then tested at z = 3 through a semi- 
analytical model of the Lya forest which adopts a lognormal Pa. 
First of all we have simulated a sample of synthetic spectra varying 



the properties of the IGM (T; To; 7) in such a way that the resulting 
mean transmitted flux matches observations. Then, we have com- 
puted Pf for each parameter set considered. Different IGM prop- 
erties affect differently the Pf, hence resulting in slightly different 
F — A relations. This provides a reconstruction method which does 
not require any assumption on the thermal and ionization state of 
the IGM. The proposed method is particularly suitable for the ex- 
traction of large scale matter density fields signals from Lya data, 
which represents the starting point for the detection of baryonic 
acoustic oscillations through quasar absortion spectra (McDonald 
& Eisenstein 2007; Slosar et al. 2009; Kitaura et al. 2010). In fact, 
this kind of studies are based on the density field reconstruction on 
scales smaller than 0. 1 Mpc along a large number of lines of sight 
more than 10 Mpc long. Such requirements are beyond the capabil- 
ity of current simulations, while can be satisfied through our fast, 
though approximate, technique. 
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